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1  Problems  Studied 

The  objective  of  the  project  is  to  extend  the  Pi’s  work  on  nonsmooth  nonlinear  equations 
from  applications  to  solvers  and  temporal  integration  to  into  nonsmooth  nonlinear  least 
squares  problems.  The  motivating  application  is  calibration  of  models  for  both  groundwater 
and  surface  water  flow.  A  secondary  objective  is  for  the  PI  to  continue  his  work  on  solvers 
and  preconditioners  for  unsaturated  flow  and  non-Darcy  models  of  saturated  flow. 

The  research  is  part  of  a  collaboration  with  a  group  at  the  Coastal  and  Hydraulics 
Laboratory  (CHL)  at  the  US  Army  Engineer  Research  and  Development  Center  (ERDC)  on 
Adaptive  Hydrology  model  (ADH),  a  production  groundwater  modeling  code,  will  be  used 
as  both  a  large-scale  testbed  for  the  algorithmic  work  in  the  project  and  as  a  source  for  new 
research  topics.  The  work  done  by  the  PI  and  his  students  has  had  a  direct  impact  on  the 
evolution  of  ADH  and  we  expect  that  the  proposed  program  of  basic  research  will  have  an 
impact  on  the  future  of  the  code.  Our  collaborators  at  ERDC  include  Charlie  Berger,  Owen 
Eslinger,  Matthew  Farthing,  Jackie  Pettway,  Stacy  Howington,  and  Chris  Kees. 

The  project  has  partially  supported  three  students.  Jill  Reese  received  her  degree  in 
2006  and  is  now  employed  at  the  Mathworks.  Corey  Winton  should  defend  his  thesis  in 
early  2010  and  formally  graduate  in  May  2010.  Winton  is  now  employed  at  the  Information 
Technology  Laboratory  at  ERDC.  Anna  Meade  is  a  new  student  and  was  supported  by  this 
project  through  May  of  2010.  she  will  be  supported  for  the  2010-2011  academic  year  directly 
by  ERDC  though  their  BAA. 
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2  Results 


We  have  published  a  paper  [4]  which  compares  several  methods  for  optimal  design  of  reme¬ 
diation  strategies.  This  effort  was  supported  by  this  project  and  the  previous  ARO  funding. 

Most  of  our  efforts  were  directed  at  reduced  order  models  for  inverse  problems  in  saturated 
flow.  Such  problems  are  a  good  first  step  to  test  ideas,  enable  us  to  integrate  the  optimization 
with  ADH,  and  may  lead  to  other  applications,  such  as  thermal  inverse  problems. 

The  computational  cost  of  an  ADH  field-scale  simulation  in  three  dimensions  is  high, 
and  it  is  important  to  reduce  the  number  of  expensive  calls  to  the  simulator.  This  is  one 
of  the  reasons  why  the  group  at  ERDC  is  not  using  methods  like  genetic  algorithms  or 
simulated  annealing  for  calibration.  Methods  based  on  nonlinear  least  squares,  however, 
miss  the  global  features  of  the  problem,  and  may  converge  to  local  minima  that  result  in 
poor  calibrations.  One  low  cost  approach  to  remedy  this  is  to  use  a  surrogate  model  or 
response  surface  to  approximate  the  expensive  objective  function,  minimize  the  surrogate 
with  a  global  optimization  method,  and  then  sample  the  expensive  function  near  one  or  more 
local  minima  of  the  surrogate. 

We  are  using  a  POD  (principal  orthogonal  decomposition)  surrogate  model  [6,8,14,16] 
which  is  a  reduced  order  model  used  in  the  fluids  control  community.  POD  uses  a  basis  which 
is  built  from  several  runs  of  the  simulator,  extracting  a  useful  low-dimensional  subspace  with 
an  eigenvalue  calculation.  We  use  a  basis  taken  from  the  sensitivity  vectors,  which  can  be 
computed  efficiently  within  finite  element  simulators  like  ADH.  The  student  involved  in  this 
part  of  the  project  is  Corey  Winton,  who  is  now  employed  at  ERDC.  Winton  has  been 
integrating  the  POD  ideas  directly  into  ADH,  and  we  have  attached  a  report  by  Winton  (see 
§  2.1)  on  the  details.  The  major  advance  in  the  last  twelve  months  is  the  integration  of  a 
well  model  into  the  POD  formulation.  Winton  should  complete  his  Ph.  D.  work  by  the  end 
of  2010.  He  has  some  significant,  but  manageable,  computing  left  to  do. 

Our  results  to  date  are  promising  [12, 13]  but  we  must  have  the  results  which  use  the 
well  model  before  publishing  in  a  journal.  The  well  model  is  required  because  the  problems 
of  interest  at  ERDC  have  integrated  flux  boundary  conditions,  and  the  well  model  maps 
those  boundary  conditions  into  pressure  boundary  conditions.  §  2.1  will  explain  the  issues 
in  detail. 

Once  the  well  model  is  in  hand,  we  have  publications  (including  Winton’s  thesis)  prepared 
and  waiting  for  it,  as  well  as  a  collaboration  with  Farthing  to  create  a  hybrid  method  that 
uses  the  POD  approximation  with  the  genetic  algorithm  from  [2]  in  the  global  phase,  and 
with  PEST  [3]  when  the  iteration  is  near  convergence.  While  genetic  algorithms  are  based 
on  heuristics,  there  are  methods  to  couple  reduced  order  models  and  heuristic  algorithms 
with  rigorous  optimization  methods  to  obtain  convergence  results  [1].  Hence,  we  expect  the 
methods  we  design  to  have  rigorous  convergence  properties. 
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2.1  Report  by  Corey  Winton 

In  the  past  year,  we  have  continued  our  application  of  Proper  Orthogonal  Decomposition 
(POD)  on  groundwater  models.  Our  goal  is  to  generate  a  reduced  model  using  POD  that 
will  allow  us  to  estimate  the  parameters  that  govern  groundwater  flow.  The  parameters 
we  are  most  concerned  with  are  hydraulic  conductivity,  the  measure  of  velocity  of  a  liquid 
through  a  porous  medium. 

This  year,  our  progress  has  focused  on  the  development  and  implementation  of  the  total 
flux  boundary  condition.  In  previous  years,  we  had  solved 


V  •  {KWh)  =  f  in  T 

(la) 

r)h 

V=0on  8YX 
on 

(lb) 

h  =  g  on  dV2 

(lc) 

We 

now  seek  to  solve  the  following  equation: 

V  •  {KWh)  =  f  in  T, 

(2) 

but  with  the  boundary  condition 

r  -  Td  u  Ttv  urQ, 

(3) 

where 

Td  =  rdl  u ...  u  TdNd, 

(4a) 

rN=  rni  u...urnjVn, 

(4b) 

Tq  =  rqi  u ...  u  r9W(J, 

(4c) 

and 

rdi  =>  h  =  adi  on  Tdi, 

(5a) 

r”,  =>  f;  = /?»,  on  r„t, 

(5b) 

/r,,!DS  =  ™rr 

(5c) 

The  solution  to  (2)  is  given  by 

h  =  h0  +  ^2jihi,  (6) 

Nq 

where  ho  and  hi  are  solutions  to  (2)  boundary  conditions  modified  as  follows: 

•  To  solve  for  h0,  we  leave  T o  and  T  tv  as  defined  in  (4),  (5),  but  force  h  =  0  on  Tq 
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dh 

•  To  compute  h,  we  set  h  —  0  on  T^,  —  =  0  on  Tjv,  and  h  —  1  on  Tq 

on 

•  The  Nq  coefficients  7 \  are  derived  from  the  following  system  of  Nq  equations: 


L 

i= 1 


r)h 

-jfo  dS  —  4>qji  j  —  l-.Nq. 


(7) 


For  the  reduced  model,  we  compute  the  solutions  hi  in  (6)  with 


Ao  +  ^  ^  A iki 


n 


h  =  fo  +  fih 


(8) 


That  is,  the  only  information  that  changes  for  each  sub-solution  hi  is  contained  in  the  vector 
f0.  The  matrix  A  is  constructed  directly  by  our  full  model  ADH.  The  vectors  /  are  derived, 
but  with  minimal  computational  effort. 

The  total  flux  formulation  with  solution  (6)  serves  several  purposes: 

•  The  total  flux  boundary  on  Tq  ensures  a  unique  solution  to  (2) 

•  We  have  more  information  in  our  basis.  Previously,  the  basis  U  was  constructed  with 


U 


,  dh 


dh 

'  dkr, 


where  k\  are  the  individual  material  conductivities  found  in  K.  Now,  we  have  consid¬ 
erably  more  information  and  can  build  the  basis 
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The  expansion  of  the  basis  is  significant  as  it  this  will  allow  the  reduced  order  model 
to  more  accurately  portray  the  physics  of  the  full  model. 


•  We  are  able  to  analytically  compute  the  gradient  for  (6)  with  minimal  additional  effort. 
We  see  from  (6)  that 


dh  _  dho  /  dhi  d^j  \ 

Skj  dkj  hV'dkj  dkj  ') 


(9) 


We  see  from  (8)  that 


dh  .  . 

— -  is  given  by 
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(10) 
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J  =  1-Nq 


We  can  find 


c>7* 

dki 


from  (7) 


#7 j 

dkk 


(11) 


Each  the  terms  in  (11)  are  straightforward  to  solve.  We  construct  a  matrix  M  that 
allows  us  to  compute  the  flux  through  any  boundary  E,  for  a  given  solution  hi.  That 


is, 

Mhi  =  [  dS. 

Jfj  on 

This  matrix  M  does  not  change  as  we  alter  the  boundary  conditions,  so  it  need  only 
be  constructed  once. 


It  should  be  noted  that  the  development  and  implementation  of  the  total  flux  boundary 
condition  is  of  interest  not  only  in  our  optimization  research,  but  also  in  the  full  ADH  code. 
This  boundary  condition  also  serves  as  a  useful  well  model  for  other  applications.  Rather 
than  present  the  well  as  a  prescribed  flux  through  a  set  of  elements,  we  can  now  solve  for 
the  accumulated  flux  through  that  boundary.  The  distinction  is  key  -  rather  than  assigning 
a  flux  at  each  element  to  simulate  a  well,  we  are  now  able  to  solve  for  the  aggregate  flux 
through  that  boundary. 


Future  Work  We  intend  to  use  the  reduced  order  model  constructed  by  POD  to  pursue 
several  optimization  goals. 

•  First,  we  intend  to  show  that  an  optimizer  using  POD  will  find  at  least  a  comparable 
local  solution  with  less  computational  effort  than  if  the  optimizer  used  the  full  model 
in  ADH. 

•  We  aim  to  use  POD  with  a  Genetic  Algorithm  (GA)  maintained  by  Matthew  Farthing 
to  find  a  global  solution  with  less  computational  effort.  A  GA  requires  a  large  number  of 
model  simulations.  We  will  demonstrate  that  the  reduced  model  constructed  by  POD 
will  allow  the  GA  to  investigate  a  large  domain  with  significantly  less  computational 
effort  than  possible  if  it  used  the  full  model. 

2.2  Other  Research  Efforts 

We  continue  to  collaborate  with  ERDC  staff  on  multilevel  solvers  (Kees)  grid  optimization 
(Eslinger),  and  maintain  discussions  on  temporal  integration  (Berger,  Howington,  Kees). 
We  have  recently  published  some  theoretical  studies  on  the  method  of  pseudo-transient 
continuation  [5,11,15,17],  which  we  have  helped  put  into  ADH.  We  have  also  done  theoretical 
work  on  nonlinear  least  squares  [7, 10] ,  which  will  have  applications  to  our  work  on  model 
calibration.  The  PI  is  also  completing  a  book  [9]  on  the  implicit  filtering  algorithm,  which 
we  have  been  using  for  design  of  remediation  systems  [4]. 
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3  Technology  Transfer  and  Collaborations  with  ERDC 

ERDC  personnel  use  the  work  of  the  Pi’s  group  in  many  ways.  Berger  has  used  the  work  on 
pseudo-transient  continuation  in  his  surface  water  work.  Our  temporal  integration  code  is 
now  completely  integrated  into  ADH,  as  is  our  work  on  preconditioning.  We  also  work  with 
Eslinger  on  optimization  methods  for  mesh  improvement. 

Over  the  term  of  the  project,  Winton  spent  8-12  weeks  at  ERDC  in  the  summers  of 
2008  and  2009.  Kelley  visited  ERDC  at  least  twice  each  year  of  the  project.  Eslinger  and 
Howington  visited  the  PI  at  NC  State  University  once  each.  The  PI  served  on  Hallberg’s 
Ph.  D.  committee. 
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